
#include <stdio.h>
#include <string.h>
#include  <stddef.h>
#include <jpeglib.h>
#include "./readbitmap.h"
#include "defs.h"
#define X265_DEPTH 8
#define ZIG(i,y,x) level[i] = dct[y][x];
#define BLOCK_COUNT 602

x265_min(int a, int b) { return a < b ? a : b; }
x265_max(int a, int b) { return a > b ? a : b; }
x265_clip3(int minVal, int maxVal, int a) { return x265_min(x265_max(minVal, a), maxVal); }

//使用的对比块来自2800与2801同一16x16块，在分块图的第17号块
const int g_t4[4][4] =
{
    { 64, 64, 64, 64 },
    { 83, 36, -36, -83 },
    { 64, -64, -64, 64 },
    { 36, -83, 83, -36 }
};

const int hdm[4][4] =
{
    { 1, 1, 1, 1 },
    { 1, 1, -1, -1 },
    { 1, -1, -1, 1 },
    { 1, -1, 1, -1 }
};

void partialButterfly4(int* src, int* dst, int shift, int line)
{
    int j;
    int E[2], O[2];
    int add = 1 << (shift - 1);

    for (j = 0; j < line; j++)
    {
        /* E and O */
        E[0] = src[0] + src[3];
        O[0] = src[0] - src[3];
        E[1] = src[1] + src[2];
        O[1] = src[1] - src[2];

        dst[0] = (int)((g_t4[0][0] * E[0] + g_t4[0][1] * E[1] + add) >> shift);
        dst[2 * line] = (int)((g_t4[2][0] * E[0] + g_t4[2][1] * E[1] + add) >> shift);
        dst[line] = (int)((g_t4[1][0] * O[0] + g_t4[1][1] * O[1] + add) >> shift);
        dst[3 * line] = (int)((g_t4[3][0] * O[0] + g_t4[3][1] * O[1] + add) >> shift);

        src += 4;
        dst++;
    }
}

//4x4DCT
void dct4_c(int **src)
{
	int i,j;
    const int shift_1st = 1 + X265_DEPTH - 8;
    const int shift_2nd = 8;

	int * coef =  calloc(4*4,sizeof(int));
	int * dst =  calloc(4*4,sizeof(int));
    int *scrData = calloc(4*4,sizeof(int));
    //int **dctData = calloc_mat_int(4,4);

    for(i = 0; i<4; i++){
		for(j = 0; j<4; j++){
			scrData[i*4+j] = src[i][j];
		}
	}

    partialButterfly4(scrData, coef, shift_1st, 4);
    partialButterfly4(coef, dst, shift_2nd, 4);

	for(i = 0; i<4; i++){
		for(j = 0; j<4; j++){
			src[i][j] = (dst[i*4+j]+16) >> 5;
		}
	}

    //return dctData;
}

void partialButterflyInverse4(int * src, int* dst, int shift, int line)
{
    int j;
    int E[2], O[2];
    int add = 1 << (shift - 1);

    for (j = 0; j < line; j++)
    {
        /* Utilizing symmetry properties to the maximum to minimize the number of multiplications */
        O[0] = g_t4[1][0] * src[line] + g_t4[3][0] * src[3 * line];
        O[1] = g_t4[1][1] * src[line] + g_t4[3][1] * src[3 * line];
        E[0] = g_t4[0][0] * src[0] + g_t4[2][0] * src[2 * line];
        E[1] = g_t4[0][1] * src[0] + g_t4[2][1] * src[2 * line];

        /* Combining even and odd terms at each hierarchy levels to calculate the final spatial domain vector */
        dst[0] = (int)(x265_clip3(-32768, 32767, (E[0] + O[0] + add) >> shift));
        dst[1] = (int)(x265_clip3(-32768, 32767, (E[1] + O[1] + add) >> shift));
        dst[2] = (int)(x265_clip3(-32768, 32767, (E[1] - O[1] + add) >> shift));
        dst[3] = (int)(x265_clip3(-32768, 32767, (E[0] - O[0] + add) >> shift));

        src++;
        dst += 4;
    }
}

void idct4_c(int ** src)
{
    int i,j;
    const int shift_1st = 2;
    const int shift_2nd = 12 - (X265_DEPTH - 8);

    int * coef =  calloc(4*4,sizeof(int));
    int * block =  calloc(4*4,sizeof(int));
    int *scrData = calloc(4*4,sizeof(int));

    for(i = 0; i<4; i++){
		for(j = 0; j<4; j++){
			scrData[i*4+j] = src[i][j];
		}
	}


    partialButterflyInverse4(scrData, coef, shift_1st, 4); // Forward DST BY FAST ALGORITHM, block input, coef output
    partialButterflyInverse4(coef, block, shift_2nd, 4); // Forward DST BY FAST ALGORITHM, coef input, coeff output

    for(i = 0; i<4; i++){
		for(j = 0; j<4; j++){
			src[i][j] = block[i*4+j];
		}
	}

}

//分成16个4x4子块的dct变换
void sub16x16_dct( int ***dct )
{
    int i,j,k;
    for ( i = 1; i < Y_Space+UV_Space; i++)
    {
        dct4_c(dct[i]);
    }
    
}

//哈达玛变换
// int *** dct4x4dc( int *** dct )
// {
//     int tmp[4][4];
//     int s01, s23;
//     int d01, d23;
//     int i;
//     int *** dctdc = calloc_mat_int3(17,4,4);
//     dctdc = dct;
//     int ** d = calloc_mat_int(4,4);
//     d = dctdc[0];

//     for( i = 0; i < 4; i++ )
//     {
//         s01 = d[i][0] + d[i][1];
//         d01 = d[i][0] - d[i][1];
//         s23 = d[i][2] + d[i][3];
//         d23 = d[i][2] - d[i][3];

//         tmp[0][i] = s01 + s23;
//         tmp[1][i] = s01 - s23;
//         tmp[2][i] = d01 - d23;
//         tmp[3][i] = d01 + d23;
//     }

//     for( i = 0; i < 4; i++ )
//     {
//         s01 = tmp[i][0] + tmp[i][1];
//         d01 = tmp[i][0] - tmp[i][1];
//         s23 = tmp[i][2] + tmp[i][3];
//         d23 = tmp[i][2] - tmp[i][3];

//         d[0][i] = ( s01 + s23 + 1 ) >> 1;
//         d[1][i] = ( s01 - s23 + 1 ) >> 1;
//         d[2][i] = ( d01 - d23 + 1 ) >> 1;
//         d[3][i] = ( d01 + d23 + 1 ) >> 1;
//     }

//     return dctdc;
// }

int *** transform1Dto3D(int * arr){
    int i,j,k;
    int *** arr3D = calloc_mat_int3(Y_Space,4,4);
    int ** arr2D = calloc_mat_int(16,16);
    for ( i = 0; i < 16; i++)
    {
        for ( j = 0; j < 16; j++)
        {
            arr2D[i][j] = arr[i*16+j%16];
        }
    }

    for ( j = 0; j < 16; j++)
    {
        for ( k = 0; k < 16; k++)
        {
            i = (j/4)*(16>>2) + k/4 + 1;
            arr3D[i][j%4][k%4] = arr2D[j][k];
        }
    }

    return arr3D;
    
}

int *** transform1Dto3D_YUV(int * arr, int * arr_uv){
    int i,j,k;
    int *** arr3D = calloc_mat_int3(Y_Space+UV_Space,4,4);
    int ** arr2D = calloc_mat_int(16,16);
    int ** arr2D_uv = calloc_mat_int(16,8);
    for ( i = 0; i < 16; i++)
    {
        for ( j = 0; j < 16; j++)
        {
            arr2D[i][j] = arr[i*16+j%16];
        }
    }

    for ( j = 0; j < 16; j++)
    {
        for ( k = 0; k < 16; k++)
        {
            i = (j/4)*(16>>2) + k/4 + 1;
            arr3D[i][j%4][k%4] = arr2D[j][k];
        }
    }

    for ( i = 0; i < 16; i++)
    {
        for ( j = 0; j < 8; j++)
        {
            arr2D_uv[i][j] = arr_uv[i*8+j];
        }
    }

    for ( j = 0; j < 16; j++)
    {
        for ( k = 0; k < 8; k++)
        {
            i = (j/4)*(8>>2) + k/4 + Y_Space;
            arr3D[i][j%4][k%4] = arr2D_uv[j][k];
        }
    }

    return arr3D;
    
}

//从文件中读取16x16宏块
int * readBlock16x16(char *path){
    FILE* fp;
    char str[100];
    int block_size=16*16;
    int *data=(int *)calloc(block_size,sizeof(int));
    fp = fopen(path, "r");
    int i=0,j=0,k=0;
    if(fp==NULL) printf("bad file path  \n");
    printf("读取宏块\n");
    while (fscanf(fp, "%s", str) != EOF)
    {
        data[i]=atoi(str);
        printf("%d\t", data[i]);
        if((i+1)%16==0){
            printf("\n");
        }
        i++;
    }
    printf("\n");
    fclose(fp);
    return data;
}

int * readBlockPosition(char *path){
    FILE* fp;
    char str[100];
    int block_size=BLOCK_COUNT*2;
    int *data=(int *)calloc(block_size,sizeof(int));
    fp = fopen(path, "r");
    int i=0,j=0,k=0;
    if(fp==NULL) printf("bad file path  \n");
    printf("读取宏块坐标\n");
    while (fscanf(fp, "%s", str) != EOF)
    {
        data[i]=atoi(str);
        // printf("%d\t", data[i]);
        // if((i+1)%2==0){
        //     printf("\n");
        // }
        i++;
    }
    printf("\n");
    fclose(fp);
    return data;
}
//从文件中读取16x16宏块
int * readBlock16x16_off(char *path,int offset){
    FILE* fp;
    char str[100];
    int block_size=16*16;
    int *data=(int *)calloc(block_size,sizeof(int));
    fp = fopen(path, "r");
    int i=0,j=0,k=0;
    if(fp==NULL) printf("bad file path  \n");
    printf("读取宏块\n");
    while (fscanf(fp, "%s", str) != EOF)
    {
        int a = rand() % offset;
        data[i]=atoi(str)+a;
        printf("%d\t", data[i]);
        if((i+1)%16==0){
            printf("\n");
        }
        i++;
    }
    printf("\n");
    fclose(fp);
    return data;
}

//打印16x16宏块
void print_16x16_block_in1D_y(int *data){
    int i=0;
    for(i=0;i<16*16;i++){
        printf("%d\t", data[i]);
        if((i+1)%16==0){
            printf("\n");
        }
    }  
}

void print_16x16_block_in1D_uv(int *data){
    int i=0;
    for(i=0;i<16*8;i++){
        printf("%d\t", data[i]);
        if((i+1)%8==0){
            printf("\n");
        }
    }  
}

void print_16x16_block_in2D(int **data, int N){
    int i,j;

    for ( i = 0; i < N; i++)
    {
        for ( j = 0; j < N; j++)
        {
            printf("%d\t", data[i][j]);
        } 
        printf("\n");
    }
    printf("\n");
    
}

void print_16x16_block_in3D(int ***data){
    int i,j;
    int ** arr2D = calloc_mat_int(16,16);
    int ** arr2D_uv = calloc_mat_int(16,8);

    for ( i = 0; i < 16; i++)
    {
        for ( j = 0; j < 16; j++)
        {
            if(i%4==0 && j%4==0){
                arr2D[i][j] = 0;
            }
            else{
                arr2D[i][j] = data[i/4*4 + j/4 + 1][i%4][j%4];
            }
        } 
    }

    int index_uv = Y_Space;
    for ( i = 0; i < 16; i++)
    {
        for ( j = 0; j < 8; j++)
        {
            index_uv = i/4*2 + j/4 + Y_Space;
            arr2D_uv[i][j] = data[index_uv][i%4][j%4];
        }
    }

    printf("DC系数：\n");
    for ( i = 0; i < 4; i++)
    {
        for ( j = 0; j < 4; j++)
        {
            printf("%d\t", data[0][i][j]);
        } 
        printf("\n");
    }
    printf("\n");
    
    printf("AC系数(0,0位置的DC系数已置0)：\n");
    for ( i = 0; i < 16; i++)
    {
        for ( j = 0; j < 16; j++)
        {
            printf("%d\t", arr2D[i][j]);
        } 
        printf("\n");
    }
    printf("\n");

    printf("UV分量：\n");
    for ( i = 0; i < 16; i++)
    {
        for ( j = 0; j < 8; j++)
        {
            printf("%d\t", arr2D_uv[i][j]);
        } 
        printf("\n");
    }
    printf("\n");
    
}

int * get_dim1_block_residual(char * path1,char *path2){
    int *block16x16_1=readBlock16x16(path1);//readBlock16x16_off(path1, 50);
    int *block16x16_2=readBlock16x16(path2);
    int *residual_data=(int *)calloc(16*16,sizeof(int));
    int i=0;
    for(i=0;i<16*16;i++){
        residual_data[i]=block16x16_1[i]-block16x16_2[i];
    }
    return residual_data;
}

int * get_dim1_block_residual_byXY(int *block16x16_1, int *block16x16_2){
    int *residual_data=(int *)calloc(16*16,sizeof(int));
    int i=0;
    for(i=0;i<16*16;i++){
        residual_data[i]=block16x16_1[i]-block16x16_2[i];
    }
    return residual_data;
}

int * get_dim1_block_residual_byXY_uv(int *block16x16_1, int *block16x16_2){
    int *residual_data=(int *)calloc(16*8,sizeof(int));
    int i=0;
    for(i=0;i<16*8;i++){
        residual_data[i]=block16x16_1[i]-block16x16_2[i];
    }
    return residual_data;
}

//获取16x16宏块的DC系数,返回包含DC系数和AC系数的17x4x4的三维数组
int *** addDCto16x16Block(int *** block){
    int i, j;
    int ** blockDC = calloc_mat_int(4,4);
    for ( i = 0; i < 4; i++)
    {
        for ( j = 0; j < 4; j++)
        {
            blockDC[i][j] = block[i*4+j+1][0][0];
        }
    }

    block[0] = blockDC;
    
    return block;
}

//量化
void quant(int ***data, int QstepDC, int QstepAC){
    int i,j,k;

    for ( i = 0; i < 4; i++)
    {
        for ( j = 0; j < 4; j++)
        {
            data[0][i][j] /= QstepDC;
        }
        
    }
    
    for ( i = 1; i < Y_Space; i++)
    {
        for ( j = 0; j < 4; j++)
        {
            for ( k = 0; k < 4; k++)
            {
                data[i][j][k] /= QstepAC;
            }
            
        }
    }

    for ( i = Y_Space; i < Y_Space+UV_Space; i++)
    {
        for ( j = 0; j < 4; j++)
        {
            for ( k = 0; k < 4; k++)
            {
                if(j==0&&k==0){
                    data[i][j][k] /= QstepDC;
                }else{
                    data[i][j][k] /= QstepAC;
                }
            } 
        }
    }
    
}

void scan_zigzag_4x4full( int *level, int **dct )
{
    ZIG( 0,0,0) ZIG( 1,0,1) ZIG( 2,1,0) ZIG( 3,2,0)
    ZIG( 4,1,1) ZIG( 5,0,2) ZIG( 6,0,3) ZIG( 7,1,2)
    ZIG( 8,2,1) ZIG( 9,3,0) ZIG(10,3,1) ZIG(11,2,2)
    ZIG(12,1,3) ZIG(13,2,3) ZIG(14,3,2) ZIG(15,3,3)
}

void scan_zigzag_4x4( int *level, int **dct )
{
                ZIG( 0,0,1) ZIG( 1,1,0) ZIG( 2,2,0)
    ZIG( 3,1,1) ZIG( 4,0,2) ZIG( 5,0,3) ZIG( 6,1,2)
    ZIG( 7,2,1) ZIG( 8,3,0) ZIG( 9,3,1) ZIG(10,2,2)
    ZIG(11,1,3) ZIG(12,2,3) ZIG(13,3,2) ZIG(14,3,3)
}

int ** getZigData(int ***data){
    int i,j;
    int ** blockZig = calloc_mat_int(Y_Space+UV_Space,16);
    
    for ( i = 0; i < Y_Space+UV_Space; i++)
    {
        if(i==0 || i>=Y_Space){
            scan_zigzag_4x4full(blockZig[i], data[i]);
        }else{
            scan_zigzag_4x4(blockZig[i], data[i]);
        }
    }
    return blockZig;
    
}

// int ** hadamard(int ** Matrix){
//     int i, j, u, v;
//     int ** HDMMatrix = calloc_mat_int(4,4);

// 	for (u = 0; u < 4; ++u) {
// 		for (v = 0; v < 4; ++v) {
//             float temp = 0;
// 			for (i = 0; i < 4; i++) {
// 				for (j = 0; j < 4; j++) {
// 					temp += (1./16)*(Matrix[i][j] * hdm[u][i] * hdm[v][j]);
// 				}               
// 			}
//             HDMMatrix[u][v] = temp;
// 		}
// 	} 
//     return HDMMatrix;
// }

void getAllresult(FILE * hFile,int filenumber1, int filenumber2, int filenumber3, int x, int y){
    char filenumberchar1[5];
    char filenumberchar2[5];
    sprintf(filenumberchar1, "%d", filenumber1);
    sprintf(filenumberchar2, "%d", filenumber2);
    char* path1 = strcat3("../dump_image/tmpfs/", filenumberchar1, ".bmp");
    char* path2 = strcat3("../dump_image/tmpfs/", filenumberchar2, ".bmp");
    int * block1 =  calloc(16*16,sizeof(int)); //存放block1的Y分量
    int * block2 =  calloc(16*16,sizeof(int)); //存放block2的Y分量
    int * block1uv =  calloc(16*8,sizeof(int)); //存放block1的UV分量
    int * block2uv =  calloc(16*8,sizeof(int)); //存放block2的UV分量
    getFullChooseArea16x16(path1, block1, block1uv, x, y, x+15, y+15);
    getFullChooseArea16x16(path2, block2, block2uv, x, y, x+15, y+15);
    int *residualy_block16x16 = get_dim1_block_residual_byXY(block1, block2);
    int *residualuv_block16x16 = get_dim1_block_residual_byXY_uv(block1uv, block2uv);
    int ***dctSub_block16x16 = transform1Dto3D_YUV(residualy_block16x16, residualuv_block16x16);

    // printf("test");
    sub16x16_dct(dctSub_block16x16);
    addDCto16x16Block(dctSub_block16x16);
    printf("两块Y分量残差：\n");
    print_16x16_block_in1D_y(residualy_block16x16);
    printf("两块UV分量残差：\n");
    print_16x16_block_in1D_uv(residualuv_block16x16);
    // 打印DC系数和AC系数
    print_16x16_block_in3D(dctSub_block16x16);
    // 对DC系数再次进行DCT变换，测试显示哈达玛聚集效果不佳
    dct4_c(dctSub_block16x16[0]);
    print_16x16_block_in3D(dctSub_block16x16);
    // 量化
    // quant(dctSub_block,7,9);
    quant(dctSub_block16x16,7,9);
    //quant(dctSub_block,100,50);+
    print_16x16_block_in3D(dctSub_block16x16);

    int ** blockZig = getZigData(dctSub_block16x16);

    char filenumberchar3[3];
    sprintf(filenumberchar3, "%d", filenumber3);
    char* path3 = strcat3("./tmpzig", filenumberchar3, "_dct");
    // save_residual_zigzig_data(path3,blockZig);
    save_allInOneFile_residual_zigzig_data(hFile, blockZig);
}



void dct_test_main(){
   
    int i,j,k;
    // int *residual_block= get_dim1_block_residual("./block.input","./block_input_1");
    // int ***dctSub_block = transform1Dto3D(residual_block);

    // sub16x16_dct(dctSub_block);
    // addDCto16x16Block(dctSub_block);

    // printf("两块残差：\n");
    // print_16x16_block_in1D(residual_block);
    // //打印DC系数和AC系数
    // print_16x16_block_in3D(dctSub_block);
    // // 对DC系数再次进行DCT变换，测试显示哈达玛聚集效果不佳
    // dct4_c(dctSub_block[0]);
    // print_16x16_block_in3D(dctSub_block);
    // //量化
    // // quant(dctSub_block,7,9);
    // quant(dctSub_block,9,13);
    // //quant(dctSub_block,100,50);
    // print_16x16_block_in3D(dctSub_block);

    // int ** blockZig = getZigData(dctSub_block);

    // save_residual_zigzig_data("./tmpzig_dct",blockZig);








    // for ( i = 0; i < 17; i++)
    // {
    //     for ( j = 0; j < 16; j++)
    //     {
    //         printf("%5d",blockZig[i][j]);
    //     }
    //     printf("\n");
        
    // }

    // int count = 0;
    // for ( i = 0; i < 17; i++)
    // {
    //     for ( j = 0; j < 16; j++)
    //     {
    //         if(blockZig[i][j]!=0){
    //             // printf("%5d\n",blockZig[i][j]);
    //             // printf("%5d\n",count);
    //             count++;
    //         }
    //     }
        
    // }
    // printf("%d\n",count);
    

    // save_residual_data("./tmp_dct",dctSub_block);

    //int num_x = {656,672,688,704,720,736,752};
    int *block16x16_postion = readBlockPosition("./blockposition1.input"); 
    int opfileflag = 0;
    FILE * hFile = fopen("./tmpall_dct","wb");

    // for ( i = 0; i < BLOCK_COUNT*2; i+=2)
    for ( i = 0; i < BLOCK_COUNT*2; i+=2)
    {
        printf("---------------------------------------------------begin---------------------------------------------------------\n");
        getAllresult(hFile, 2282,2281,opfileflag,block16x16_postion[i],block16x16_postion[i+1]);
        printf("---------------------------------------------------end---------------------------------------------------------\n");
        opfileflag++;
    }
    fclose(hFile);
    

    // for ( i = 0; i < 894; i++)
    // {
    //     printf("%5d",block16x16_postion[i]);
    // }
    
    // getAllresult(2277,2276,0,656,208);

    // int filenumber1 = 2272;
    // int filenumber2 = 2273;
    // char filenumberchar1[5];
    // char filenumberchar2[5];
    // sprintf(filenumberchar1, "%d", filenumber1);
    // sprintf(filenumberchar2, "%d", filenumber2);
    // char* path1 = strcat3("../dump_image/tmpfs/", filenumberchar1, ".bmp");
    // char* path2 = strcat3("../dump_image/tmpfs/", filenumberchar2, ".bmp");
    // int * block1 =  calloc(16*16,sizeof(int));
    // int * block2 =  calloc(16*16,sizeof(int));
    // getFullChooseArea16x16(path1, block1, 576, 263, 591, 278);
    // getFullChooseArea16x16(path2, block2, 576, 263, 591, 278);
    // int *residual_block16x16 = get_dim1_block_residual_byXY(block1, block2);
    // int ***dctSub_block16x16 = transform1Dto3D(residual_block16x16);

    // // printf("test");
    // sub16x16_dct(dctSub_block16x16);
    // addDCto16x16Block(dctSub_block16x16);
    // printf("两块残差：\n");
    // print_16x16_block_in1D(residual_block16x16);
    // //打印DC系数和AC系数
    // print_16x16_block_in3D(dctSub_block16x16);
    // // 对DC系数再次进行DCT变换，测试显示哈达玛聚集效果不佳
    // dct4_c(dctSub_block16x16[0]);
    // print_16x16_block_in3D(dctSub_block16x16);
    // //量化
    // // quant(dctSub_block,7,9);
    // quant(dctSub_block16x16,9,13);
    // //quant(dctSub_block,100,50);
    // print_16x16_block_in3D(dctSub_block16x16);

    // int ** blockZig = getZigData(dctSub_block16x16);

    // int filenumber3 = 0;
    // char filenumberchar3[3];
    // sprintf(filenumberchar3, "%d", filenumber3);
    // char* path3 = strcat3("./tmpzig", filenumberchar3, "_dct");
    // save_residual_zigzig_data(path3,blockZig);
}
